/**
 * 
 */
package ua.lviv.franko.mlsFunctions.threeDimension;

import java.util.ArrayList;

import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.threeDimension.weightFunctions.CubeSupport;
import ua.lviv.franko.mlsFunctions.threeDimension.weightFunctions.IWeightFunction3D;
import ua.lviv.franko.visualisation.Visualise3D;

/**
 * @author marko
 * 
 */
public class BaseFunction {

	protected int					PolinomialDegree;
	protected Particle				xi;
	protected IWeightFunction3D		func;
	protected HashMatrixLUStyle			A;
	protected double[]				b;
	protected ArrayList<Particle>	NodeList;
	protected double				EPS	= 0.0000001;

	public BaseFunction(int PolinomialDegree, Particle xi, IWeightFunction3D func) {
		this.PolinomialDegree = PolinomialDegree;
		this.xi = xi;
		this.func = func;
	}

	public double calculate(double x, double y, double z) {
		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}
		if ((xi.zi - xi.rho > z) || (xi.zi + xi.rho < z)) {
			return 0;
		}
		ArrayList<Double> BaseVector = getBaseVector(x, y, z, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new HashMatrixLUStyle(RowType.ARRAY_TYPE, n);
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, z, particle);
		}

		double val = func.calculate(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculate(x,
																			// xi.xi,
																			// xi.rho);
		BaseVector = getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			LUFactorization lu = new LUFactorization(A);
			double[] res = lu.calculate(b);

			BaseVector = getBaseVector(xi.xi, xi.yi, xi.zi, xi);
			for (int i = 0; i < b.length; i++) {
				returnValue += BaseVector.get(i) * res[i];
			}
		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
	}

	public double calculateDerivatyX(double x, double y, double z) {
		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}
		if ((xi.zi - xi.rho > z) || (xi.zi + xi.rho < z)) {
			return 0;
		}
		ArrayList<Double> BaseVector = getBaseVector(x, y, z, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new HashMatrixLUStyle(RowType.ARRAY_TYPE, n);
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, z, particle);
		}

		double val = func.calculate(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculate(x,
																			// xi.xi,
																			// xi.rho);
		BaseVector = getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			LUFactorization lu = new LUFactorization(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyX(x, y, z, xj.xi, xj.yi, xj.zi, xj.rho);// func.calculateDerivaty(x,
																					// xj.xi,
																					// xj.rho);//
																					// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj.zi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// альфа на вектор Р і
														// вагу ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// число на вектор Р
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyX(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculateDerivaty(x,
																				// xi.xi,
																				// xi.rho);
																				// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi.zi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			if (key) {
				double[] res2 = lu.calculate(b);
				BaseVector = this.getBaseVectorDerivatyX(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res[j];
				}
				BaseVector = this.getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// return 0.0;

	}

	public double calculateDerivatyY(double x, double y, double z) {
		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}
		if ((xi.zi - xi.rho > z) || (xi.zi + xi.rho < z)) {
			return 0;
		}
		ArrayList<Double> BaseVector = getBaseVector(x, y, z, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new HashMatrixLUStyle(RowType.ARRAY_TYPE, n);
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, z, particle);
		}

		double val = func.calculate(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculate(x,
																			// xi.xi,
																			// xi.rho);
		BaseVector = getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			LUFactorization lu = new LUFactorization(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyY(x, y, z, xj.xi, xj.yi, xj.zi, xj.rho);// func.calculateDerivaty(x,
																					// xj.xi,
																					// xj.rho);//
																					// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj.zi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// альфа на вектор Р і
														// вагу ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// число на вектор Р
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyY(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculateDerivaty(x,
																				// xi.xi,
																				// xi.rho);
																				// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi.zi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			if (key) {
				double[] res2 = lu.calculate(b);
				BaseVector = this.getBaseVectorDerivatyY(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res[j];
				}
				BaseVector = this.getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// return 0.0;

	}

	public double calculateDerivatyZ(double x, double y, double z) {
		if ((xi.xi - xi.rho > x) || (xi.xi + xi.rho < x)) {
			return 0;
		}
		if ((xi.yi - xi.rho > y) || (xi.yi + xi.rho < y)) {
			return 0;
		}
		if ((xi.zi - xi.rho > z) || (xi.zi + xi.rho < z)) {
			return 0;
		}
		ArrayList<Double> BaseVector = getBaseVector(x, y, z, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new HashMatrixLUStyle(RowType.ARRAY_TYPE, n);
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, z, particle);
		}

		double val = func.calculate(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculate(x,
																			// xi.xi,
																			// xi.rho);
		BaseVector = getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			LUFactorization lu = new LUFactorization(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyZ(x, y, z, xj.xi, xj.yi, xj.zi, xj.rho);// func.calculateDerivaty(x,
																					// xj.xi,
																					// xj.rho);//
																					// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj.zi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// альфа на вектор Р і
														// вагу ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// число на вектор Р
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyZ(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculateDerivaty(x,
																				// xi.xi,
																				// xi.rho);
																				// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi.zi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			if (key) {
				double[] res2 = lu.calculate(b);
				BaseVector = this.getBaseVectorDerivatyZ(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res[j];
				}
				BaseVector = this.getBaseVector(x, y, z, new Particle(x, y, z, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}
		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// return 0.0;

	}

	public ArrayList<Double> getBaseVector(double x, double y, double z, Particle xi) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; i + j <= this.PolinomialDegree; j++) {
				for (int k = 0; i + j + k <= this.PolinomialDegree; k++) {
					L.add(Math.pow((xi.xi), i) * Math.pow((xi.yi), j) * Math.pow(xi.zi, k));
				}
			}
		}

		return L;
	}

	public ArrayList<Double> getBaseVectorDerivatyX(double x, double y, double z, Particle xi) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; i + j - 1 <= this.PolinomialDegree; j++) {
				for (int k = 0; i + j + k <= this.PolinomialDegree; k++) {
					L.add(i * myPow((xi.xi), i - 1) * myPow((xi.yi), j) * myPow(xi.zi, k));
				}
			}
		}
		return L;
	}

	public ArrayList<Double> getBaseVectorDerivatyY(double x, double y, double z, Particle xi) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; i + j <= this.PolinomialDegree; j++) {
				for (int k = 0; i + j + k <= this.PolinomialDegree; k++) {
					L.add(j * myPow((xi.xi), i) * myPow((xi.yi), j - 1) * myPow(xi.zi, k));
				}
			}
		}
		return L;
	}

	public ArrayList<Double> getBaseVectorDerivatyZ(double x, double y, double z, Particle xi) {
		ArrayList<Double> L = new ArrayList<Double>();
		for (int i = 0; i <= this.PolinomialDegree; i++) {
			for (int j = 0; i + j <= this.PolinomialDegree; j++) {
				for (int k = 0; i + j + k <= this.PolinomialDegree; k++) {
					L.add(k * myPow((xi.xi), i) * myPow((xi.yi), j) * myPow(xi.zi, k - 1));
				}
			}
		}
		return L;
	}

	public void setNodeList(ArrayList<Particle> nodeList) {
		NodeList = nodeList;
	}

	private void calculateAlfa(double x, double y, double z, Particle xi) {
		ArrayList<Double> BaseVector = getBaseVector(x, y, z, xi);// getBaseVector(x,
																	// xi);
		int n = BaseVector.size();
		for (int i = 0; i < n; i++) {
			double val = BaseVector.get(i) * func.calculate(x, y, z, xi.xi, xi.yi, xi.zi, xi.rho);// func.calculate(x,
																									// y,
																									// xi.xi,
																									// xi.yi,
																									// xi.rho);//func.calculate(x,
																									// xi.xi,
																									// xi.rho);
			for (int j = 0; j < n; j++) {
				A.setValue(i, j, A.getValue(i, j) + val * BaseVector.get(j));
			}
		}
	}

	private double myPow(double x, double pow) {
		if (pow < 0) {
			return 0.0;
		} else {
			return Math.pow(x, pow);
		}
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		ArrayList<Particle> nodes = new ArrayList<Particle>();
		// параметри області:
		double ax = 0, bx = 1, ay = 0, by = 1, az = 0, bz = 1;

		int n = 15, m = 15, mn = 15;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m, hz = (bz - az) / (double) mn;
		;

		// вузли+функції:
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				for (int k = 0; k < mn; k++) {
					Particle p = new Particle(ax + i * hx, ay + j * hy, az + k * hz, 0.5);
					nodes.add(p);
				}
			}
		}
		final BaseFunction ff = new BaseFunction(1, new Particle(0.5, 0.5, 0.5, 0.5), new CubeSupport(new CubicSpline()));
		ff.setNodeList(nodes);

		Visualise3D.show(new Mapper() {

			@Override
			public double f(double x, double y) {
				return ff.calculate(x, y, 0.6);// (x, y);
			}
		}, new Range(0.0, 1.0), new Range(0.0, 1.0));

	}

}
